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i-Q " Abstract 
^. 

pLn . The characterization of particle diffusion is a classical problem in physics and probability 

■ theory. The field of microrheology is based on experiments in which microscopic tracer beads 
P\| . are placed into a non-Newtonian fluid and tracked using high speed video capture. The 

modeling of the behavior of these beads is now an active scientific area which demands multiple 
stochastic and statistical methods. 

We propose an approximate wavelet-based simulation technique for two classes of con- 
tinuous time anomalous diffusion models, the fractional Ornstein-Uhlenbeck process and the 
fractional generalized Langevin equation. The proposed algorithm is an iterative method that 
provides approximate discretizations that converge quickly and in an appropriate sense to the 
^ , ' continuous time target process. As compared to previous works, it covers cases where the 

natural discretization of the target process does not have closed form in the time domain. 
Moreover, we propose smoothing procedures as to speed the time domain decay of the filters. 

>' 
l> 

^ ; 1 Introduction 

. The characterization of particle diffusion is a classical problem in physics and probability theory 

■ dating back to the early work of Einstein on Brownian motion. The case of Newtonian fluids, 
CN ■ which includes water, is now well-understood; however, the case of complex fluids is an active 

area of experimental and theoretical inquiry. Of particular interest are fluids of biological origin 
such as mucus, which is both highly viscous due to the presence of compounds such as mucin, a 
^ ■ high molecular weight protein that enters into the composition of human mucus, and elastic due 
^ . to the cross-linking of the mucin. Models of diffusion in biological fluids have been developed for 
many pharmaceutical and medical applications (e.g., Lai et al. (2009), Suk and Lai (2009), Suh 
et al. (2005), Dixit et al. (2008)). The related field of microrheology is based on experiments in 
which tracer beads, whose radii range from tens of nanometers to micrometers, are placed into the 
complex fiuid and tracked at millisecond sampling rates using high speed video capture and light 
microscopy (e.g.. Mason and Weitz (1995)). Tens of those beads can be tracked at one time, and 
the data is preprocessed using particle tracking software to yield the position X{t), as a function 
of time, of the individual diffusing particles. The modeling of the behavior of these beads calls for 
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a wide range of stochastic and statistical methods, including simulation (see Didier et al. (2012) 

and references therein). 

In the physics literature, a position process is called diffusive if its second moment satisfies 
Einstein's Brownian diffusion law, i.e., 

EX^{t) r^e, oo, (1.1) 

with a = 1. However, if a ^ 1, then it is said to be anomalously diffusive. The subcases 
of < a < 1 and a > 1 are named sub- or super diffusive, respectively (see Kou (2008), and 
references therein). The primary focus of this paper is to provide a fast, approximate wavelet- 
based method for a certain class of anomalous diffusion models. We approach this problem with 
wavelet methods similar to those found in the papers of Didier and Pipiras (2008, 2010), which 
in turn are related to many other works including Meyer et al. (1999), Zhang and Walter (1994), 
Sellan (1995), Pipiras (2004). 

For water, the velocity V of an embedded free particle satisfies the Langevin equation. The 
stationary solution for V gives the well-known Ornstein-Uhlenbeck process (OU). In order to 
encompass the broad class of non-Newtonian, viscoelastic fluids such as many biological fluids, 
one cannot assume that the influences of the various fluid molecule impacts are independent of 
recent activity by the free particle. This non-Markovian property is modeled by a memory kernel, 
which we denote by r(t). The generalized Langevin equation (GLE) for the velocity process of a 
free particle is 

m^F(t) = -Cy" T{t-s)V{s)ds + F{t), (1.2) 

where m is the particle mass, ( is the friction constant, ks is the Boltzmann constant, and F(t) 
is a stationary Gaussian process with autocorrelation function 

EF{t)F{s) = kBTCT{\t-s\), (1.3) 

where r is the temperature. The special case of the Langevin equation is obtained by setting F 
to a Dirac delta distribution. In this paper, wc are particularly interested in the (subdiffusive) 
fractional GLE (fGLE), since its correlation structure is now well-studied (see Kou (2008)). This 
corresponds to setting F{t)dt to dBffit) up to a constant, where BH{t) is a fractional Brownian 
motion (fBm). The latter is the unique Gaussian, if-self-similar, stationary increment process 
(see, for instance, Taqqu (2003)). In the framework of the fGLE, the memory kernel satisfies 

V{t) = 2H{2H -l)\t\^^-'^, tj^O, ^<H<1. (1.4) 

The fractional Ornstein-Uhlenbeck process (fOU) is another model for anomalous diffusion of 
interest in this paper and is the a.s. continuous solution to the fBm-driven Langevin equation 

dV{t) = -CV{t)dt + adBH{t), t>0, 0<H <1 (1.5) 

(see Cheridito et al. (2003), Prakasa Rao (2010)). Its simulation is interesting in its own right, 
since the fOU is a model for both sub- and superdiffusion. In addition, simulation of the fOU can 
be viewed as a step towards the simulation of the full fGLE, both for analytical convenience and 
due to their similar correlation structures. To view qualitative differences between the processes, 
sample paths of the OU, fOU, and fGLE processes are shown in Figure 1. 

Our proposed simulation procedure is based on a wavelet analysis of the velocity process V. 
This analysis encompasses three components: 
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(i) a wavelet-based decomposition of V{t); 



(ii) a sequence of stationary discrete time processes Vj = {Vj^n}nez, 

oo 

yj,k= E 9j,n^k-n, &}'~ WN(0,1), (1.6) 



that can be thought of as discrete approximations of the continuous time process V at scale 

(iii) a Fast Wavelet Transform (FWT)-like algorithm relating Vj across different scales, as well 
as a fast convergence of 2^/'^Vj to the continuous time process F as j ^ oo (small scales). 

The simulation procedure itself can be directly expressed as components (ii) and (Hi). More 
details about the wavelet-based decomposition (i), the theoretical backbone of the simulation 
technique, can be found in Didier and Pipiras (2008). We now explain the components (ii) and 
(Hi) in more detail including what needs to be calculated to implement the simulation algorithm. 

The simulation method relies on the Cramer- Wold Fourier domain representation of the sta- 
tionary velocity process 

Vit) = [ e''-g{x)B{dx), (1.7) 

JR 

where 'g{x) G L^(R) is called a spectral filter and B(dx) is a complex valued Brownian measure 
dominated by the Lebesgue measure. Throughout the paper, the Fourier transform of either 
a discrete or continuous time function/filter / is denoted by /. In either the discrete (1.6) or 
continuous (1.7) time setting, the autocovariance function of the process in question is given by 
the inverse Fourier transform of the spectral density, i.e., of the functions |^j(a;)p and [^(a;)^, 
respectively. The simulation method relies on designing an appropriate sequence of discrete time 
filters gj, j = 0,1,..., J, where J is the finest approximation scale chosen. The filters are then 
used to recursively generate an induced sequence of discrete approximations Vj, at scale j G N. 
The choice of the sequence of discretization filters gj is the key requirement. Heuristically, 

Gj-(2-^x) := f« 1, i^oo, xeR, (1.8) 

where gj{x) is periodically extended to M. Intuitively, the discretization filter approximates, up 
to a scaling factor, the continuous time filter as the scale becomes finer and finer. The simulation 
procedure can then be briefly described as the following FWT-like algorithm. 

{Wi): at scale/step 0, generate via an exact method (e.g., Circulant Matrix Embedding) one first 
approximation sequence Vq; 

(W2): at scale/step j G N, and given an approximation sequence Vj, obtain the next discrete 
approximation T^+i at scale/step j + 1 via the relation 

Vj+i = Uj* t2 Vj + Vj* t2 £j, 

where {sj} is a noise sequence, and (t2 ^Ofc = ^■,fe/2l{cvcn k} + 01{odd k} is the upsampling by 
factor 2 operator. The wavelet filters Uj and Vj are defined in the Fourier domain as 

%(^) := %7rT^^^)' Vj{x) ■-gj+i{x)v{x), (1.9) 
9j j 
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where u, v are the conjugate mirror filters (CMF) of an underlying wavelet multiresolution analysis 
(MRA). For the theoretical purposes of this paper, we use a Meyer MRA, due to the compact 
Fourier domain support of both the scaling and wavelet functions (see Mallat (1999), chapter 7). 

The algorithm works because at step j + 1, the FWT annihilates the correlation structure gj 
of the approximation at scale j (see expression (1.6)) and replaces it with a new, pre-chosen 

correlation structure gj-\-i- The properties of the CMFs u and v play an important role, which we 
can explain heuristically. By taking Fourier transforms on both sides of (1.6), Vj{x) = gj{x)^{x), 
where ^ is a white noise sequence and ^(x) is its Fourier transform. Then 

^^■+1^^) = T^^(^)Vj{2x)+gj+,{x)v{x)e{2x), 

where e is another independent white noise sequence. Therefore, Vjj^i{x) = gjj^i{x){u{x)^{2x) + 
v{x)£{2x)). Since u and v are CMFs, then the term u[x)^{2x) +v{x)e{2x) is itself distributed as 
white noise. Thus, in law, Vjj^i{x) = gjj^i{x)r}{x), where rj is another white noise sequence; see 
also Pipiras (2005) for details. Moreover, we can show that 

2^/'K7,L2^tJ « J^^. (1-10) 

in a suitable sense, where \_x\ denotes the integer part of x G M (see Proposition A.l). If (1.10) 
holds, then the discrete sequences {Vj',n} can be used as a discrete approximation to V ^ which 
in turn leads to a natural approximation of the position process X(t) as a Riemann sum (see 
Proposition A. 2). Why should the property (1.10) hold? An approximation term can be expressed 
as the integral 

y^ „= f V{t)^j{t-2-^n)dt, 

Jr 

where $j can be suitably defined in the Fourier domain as 

8j(x) := Gj{2-jx)2-^/^'${2-^x), (1.11) 

and (j) is the scaling function of the underlying MRA. Since 2~"^/^0(2~'^x) 2~'^/^ as J ^ oo, 
then 8j(x) 2--^/2^(2--^x). Therefore, 

2-^/Vj„ = ^ / F(x)e-^2"^"$j(x)dx ~ 7^ / F(x)e-^^2-^n^^ ^ V{2-^n), 
27r 7]R 2TT 7k 

as desired. 

The choice of the function Gj{x) plays a key role in the implementation of the simulation 
algorithm. The case of the OU process is quite special; its discretization based on time intervals 
of any size leads to an AR(1) sequence. In the simulation framework of this paper, the specification 

Gj(x) = 2-^^—— . \ / (1.12) 

' (C + 2Jzx)-i ^ ' 



is a natural choice for the OU process, since (1 — e~^^ ''e~'^)~^ is a spectral filter of the associated 
discrete time AR(1) sequence. For the sake of exposition, note that the heuristic relation (1.8) 
is indeed satisfied in a pointwisc sense. Some exact and analytically convenient discretization 
schemes such as (1-12) are discussed in Didier and Pipiras (2008). However, simple discretization 
schemes are not available for most classes of stochastic processes. As a rule, the discretization of a 
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continuous time process leads to a substantially more intricate expression for the spectral density, 
such as for fBm and the fGLE. Though expressions for the covariance structure are available for 
the fGLE, they do not appear in closed form and will require numerical methods. 

The main contribution of this paper is to go beyond Didier and Pipiras (2008) in the sense 
of proposing a simple method to obtain discretization spectral filters (i.e., Gj) for the velocity 
process of the fGLE and the fOU. Unlike the OU process, the spectral densities of the associated 
discrete time processes do not have closed form Fourier series representations and exhibit fractional 
behavior at the origin. The method involves truncation in the Fourier domain of the continuous 
time spectral filter, and generates non-causal filters. We show that the proposed filters have 
quadratic decay, which is efficient for computational approximation. Moreover, we also propose 
a technique to smooth the filters in the Fourier domain to accelerate their time domain decay. 
Thus, the discretization sequences {V^-,n} obtained are not exact discretizations (at scale j) of 
the continuous time process V, but we show that property (1.10) still holds. The simulation 
procedure shares the positive properties of Didier and Pipiras (2008), such as: it is computationally 
fast, potentially reaching complexity 0{N), since it is based on a Fast Wavelet Transform-like 
algorithm; it provides approximations that converge uniformly over compact intervals almost 
surely; the convergence speed is exponentially fast and depends on the sample path smoothness 
of the limiting process. It is also iterative both intensively and extensively. In other words, a 
generated approximation at scale J over a compact interval [0, T], T G R+, can be used to generate 
a finer approximation at scale J+1 over [0, T] or some expanded interval [0, T + x] , X ^ ■ Our 
simulation procedure also serves to simulate the position process X with the same convergence 
rate as the velocity process V over compact intervals. Moreover, the method is not intrinsically 
Gaussian, although Gaussian processes are the primary focus of this paper. 

The paper is structured as follows. In Section 2, we introduce the Fourier domain integral 
representations for the velocity processes V for the fGLE and the fOU process, then we develop 
the discretization filters that enter into the simulation procedure. In Section 3, we evaluate 
the accuracy of the wavelet-based simulation method in comparison to other, exact methods. 
Appendices A and B contain all the proofs and auxiliary results, respectively. Appendix C contains 
the tables with the simulation results, while Appendix D shows a study of the numerical accuracy 
of the computational techniques. For clarity, pseudocode for the simulation method is provided 
in Appendix E. 

2 Filters in the spectral domain 

In this section, we propose approximate discretization filters with the purpose of simulation. The 
focus is on the fGLE and fOU, but the OU process will be revisited frequently in order to contrast 
exact and approximate discretization procedures. 

2.1 Fourier domain representations 

First, we express the integral representations which will be the basis for the construction of the 

simulation procedures for each class of processes. It is convenient to rewrite the velocity process 
y as a Fourier domain stochastic integral with respect to a Brownian measure. So, we define Bi, 
B2 as two real- valued Brownian motions and B{dx) = Bi{dx) + iB2{dx) as the induced random 
measure satisfying 

B{-dx) = -B{dx) a.s., E\B{dx)\'^ = dx. (2.1) 
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Remark 2.1 In this paper, wc use the parameter S to represent the fractional behavior of a 
spectral density around the origin. In other words, 

\g{x)\'^x-'\ -1<S<1. (2.2) 

If (5 > or (5 < , then the process is said to be long range dependent or antipersistent, respectively. 
As shown in this section, both the fOU and the fGLE have spectral densities that satisfy (2.2). 
Moreover, in both cases we will define 

d = H -1/2, (2.3) 

where H is the Hurst parameter of the driving fBm. Note, however, that for the fOU, 6 = d, 
whereas for the fGLE, 6 = —d. In both cases, a in (1.1) and 5 are connected by means of the 
relation a = 1 + 2(5 (see Didier et al. (2012)). 

The fOU process admits the spectral representation 



^__/r(2d + 2)sin(7r(d + l/2)) /■ 1 



V{t) = ad^ ' / e^*^ , \x\-''B(dx), --<d<-, (2.4) 

from Proposition 2.5 of Cheridito et al. (2003), pp. 5-8. 

The integral representation of the OU process can be obtained by setting d = in (2.4), where 
(2.3) holds and H is as in (1.5). The velocity process for the fGLE (1.2), (1.4) can be represented 
as 

^ / ,\ / Jtx J- l^lrfr^/ ' ^ . , , 1 



Vit) = c{d) / e'^-—^-——-——\x\''B{dx), < d < -, (2.5) 

where j5 = 1 + 2d, and the constants are defined by 

7o = a^ + 62^ 71 = 26m, 7 = c{d) = ^J2C,kBT V{2d + 2) sin(7r(d + 1/2)) 

and a = CX{2d + 2) sin(7r(d + 1/2)), h = Cj:{2d + 2)cos(7r(d + 1/2)). Note that C > 0, and 
< d < 1/2, so a > and 6 < 0. 

The expression (2.5) can be established by following the techniques of Kou (2008), Theorem 
2.1. One starts by expressing the velocity process of a fBm-driven free particle (1.2), (1.4) in 
terms of a pathwise defined Riemann-Stieltjes integral 

V{t) = ^/2CkBT f r{t - u)BH{du). (2.6) 

The time domain filter r is given through the inverse Fourier transform 

r{t) = ^ f ~, ^ e-^*^dx, 

271" Jr CKJj (x) — imx 

where K^{x) = \x\^~'^^T[2H + l)(sin(if7r) — isign(a;) cos(i?7r)). To obtain (2.5), rewrite the 
integrand (2.6) as a fractional integral as in Pipiras and Taqqu (2000) or, equivalently, adapt the 
expression for the spectral density of V developed in Kou (2008), p. 524. We arrive at a Fourier 
domain integral with respect to a measure (2.1). The spectral filter is 

(^k(x) — imsign(x)|x|^+^'^ ' 

up to a constant; this leads to the filter in (2.5) by elimination of the imaginary part. 

Expression (2.4) shows that the fractional parametrization 5 (see (2.2)) of the fOU encompasses 
the full range in (2.2), whereas the fGLE is necessarily antipersistent. For the sake of illustration, 
the correlation structure of the fGLE is depicted in Figure 2. Its autocovariance function displays 
the characteristic fluctuations associated with antipersistence. 
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Figure 2: fGLE, correlation structure up to a multiplicative constant (d = 0.10,0.25,0.35,0.45, 
= 2, m = 1). Left plot: spectral density. Right plot: autocovariance function. 



2.2 Wavelet decomposition and the choice of a discretization filter gj 

In this subsection, we will present the wavelet-based simulation method and summarize previous 
results relevant to the current situation. In addition, we will present the computational details of 
our proposed implementation for the simulation procedure. Before proceeding to details, we will 
present sufficient conditions for the underlying wavelet decomposition to hold and give intuitive 
explanations for these assumptions. 

In order to construct a wavelet-based decomposition of the velocity process V, the following 
technical assumptions must be met (see Didier and Pipiras (2008)). 



Assumption 1: 
Assumption 2: for any j G Z, 
Assumption 3: for any jo G Z, 



loc 



loci 



max max sup sup 

p=-l,lfe=0,l,2j>jo \x\<iTr/3 



d''{Gj{x))P 



< oo. 



(2.7) 



(2.9) 



Assumption 4: for large \x\ 



Assumption 5: for large j, 



d^'g{x) 



dx^ 



const 

<r-^, fc = 0,l,2. 



\Gj{0) - 1| < const 2"^. 



(2.10) 



(2.11) 



Assumption 1 pertains only to g and is clearly satisfied by (2.4), (2.5). In order to clarify and 
explicitly incorporate processes with fractional spectral densities, Assumption 4 is replaced in 
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Appendix A by the slightly modified Assumption 4 . On the other hand, Assumptions 2, 3 and 5 
pertain Gj, the choice of the discretization filter gj. As explained in Section 1, the discretization of 
continuous time processes will in general generate discrete time processes whose spectral densities 
have intricate analytical forms. So, a natural question is whether one could construct converging 
discretizations by means of a simple method that applies to a wider class of stochastic processes, 
in particular, the fGLE. Indeed, this can be done by developing non-causal discretization filters 
gj in three elementary steps: 

(t.l) extend the truncated function 5(2:)l[_7r,7r) periodically to M implying that 'gj stems directly 

from g; 

(t.2) modify the resulting function with rescaling terms (e.g., 2^) so that relation (1.8) holds; 

{t.3) smooth the resulting function at — tt and tt to speed up the time domain decay of the filter 
in theory and computational practice. 

Step (t.l) replaces exact discretization filters such as those for the AR(1). Step (t.2) is necessary to 
ensure the coherence between the discretization and the limiting process. Step (t.3) minimizes the 
border effect, the consequence of the truncation of infinite-length filters which is always present 
in computational practice. The method described in steps (t.l) — (t.3) produces a sequence of 
discrete time processes Vj = gj * £ which is approximate, since gj is picked for analytical and 
computational convenience. 

We first look at the truncated procedure and filters obtained via steps (t.l) and (t.2). The 
idea of constructing gj{x) by truncating g{x) at ibvr and extending the function periodically to 
M has the obvious advantage of being a simple method for obtaining a discretization. Even 
though it will typically create filters with discontinuities at ±7r in the derivatives of the function 
gj{x), it is still possible to show that the filters exhibit good theoretical decay under assumptions. 
Furthermore, one can show that these approximate discretizations still converge exponentially 
fast to the correct limiting process. Mathematically, this requires replacing Assumption 3 with a 
set of weaker conditions. The details are described and established in Appendix A. 



Table 1: Spectral density components and associated truncated filters 
process IffO^IF dj (•^) 

ou, fou \e + x^\-^ (c^ + 22^x2);i/' 

fOU, fGLE Ixp'^ 2^'^\x\^ 

fGLE |7o + 7ik|^ + 72|a;|^^|-i {lo + li2^^\xf + 722^^1^ \x\^f^)p 



We now propose truncated filters for the processes of interest. Table 1 contains the truncated 
filters generated based on components of the spectral filters for the OU, fOU and fGLE processes. 
The subscript p denotes periodic extension beyond the domain [— 7r,7r). In all cases, we deal with 
purely real filters for analytical simplicity. The associated functions (1.8) are 

cu) (C^ + 2^^-.x^);^/- 

9jd\rr\d \^\d 

GjA^) = ^ = (2-13) 

r . ._ (70 + 7i2^-^|x|/^ + 722^^-/^|xp/^);^/^ 
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frequency time 

Figure 3: OU process, filters at j = 1 ((^ = 1, no zero moments). Left: exact versus non-smoothed 
liigh-pass filter, spectral domain. Right: non-smoothed high-pass filter, time domain. 



for x G M. 

Example 2.1 In order to contrast the exact and approximate discretization filters we revisit the 
well-known case of the OU process. The rationale behind the proposed gj{x) is simple. Multiplying 
the argument of a function by a number less than 1 has the geometric effect of "stretching" 
the original function away from the origin. Therefore, we have that gj{2~^x) = (C^ + 
(pointwise) for any fixed x G M and large enough j, as desired. In other words, the discrete time 
process whose spectral density is gj{x) is, indeed, an approximation to the OU process in the 
sense explained in Section 1. The difference between the truncated filter and the exact filters 
(1.12) can be seen in the Fourier domain in Figure 3, left plot. As expected, even though gj{x) is 
continuous at ±7r, the same is not true for its first derivative. The left derivative of gj at tt is 



^J-^ ^ ^^2 _^22i7r2)3/2■ 
Moreover, since g^{x) = l/\/C'^ + 2^^ (27r — a;)^ for x G [tt, 27r), then the right derivative of gj at 
vr is gj _,_(7r) = —gj _{'k) by periodic extension. 



Although the theoretical guarantees on the asymptotic decay of the truncated filters are good, 
the discontinuities at ±tt in the derivatives of the function gj{x) will create ripples in the time 
domain expression of the associated filters. This can be seen for the OU process filters in Figure 
3, right plot. In computational practice, one could ask whether it is possible to eliminate or at 
least to minimize the ripples one observes in the plots. Step (t.3) consists of addressing this issue 
by smoothing the kinks of gj{x) at {{2k + l)vr}/;gg. There is clearly more than one way to do 
this. Table 2 contains the proposed smoothed filters. Denote the discretization filters in the first 
through third rows of Table 1 by gj^(^{x), gj^dix), gj,'y,dix), respectively. Thus, the proposed filters 
for the OU, fOU and fGLE processes are gj^(^{x), gj^c{x)gj^d{x) and gj^Q{x)gj^^^d{x)^ respectively. 
In contrast with exact discretization filters (see Figure 4), they generally give rise to non-causal 
low- and high-pass filters Uj and Vj, which is inconsequential for simulation purposes. 
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Table 2: Spectral density components and associated truncated-smoothed filters 



process \9{x)\ 9j{^) > 

- ^, where v > ^, - ^ ^ 



OU,fOU \e + xr' .rSi-tfiv. V^-&,where.>0,^>^ 



„„ „ I- £ign(£) £*(i)£ \ 2 

fou,fGLE '^^s .";:t7^ ^ -vi^ 



Ztt 

(70+712^" kl''+7222i/5|a;|2/3)^/ 








exact, AR{1} 

approximate 






H 





-0.3 

-40 -30 -20 -10 




time 



10 20 30 40 



Figure 4: OU process, exact versus approximate high-pass filters at j = 1 (C = 1, 4 zero moments). 
Left: spectral domain. Right: time domain. 



Example 2.2 Multiplying the original truncated discretization filter for the OU process by a 

term of the form creates in the former two global minima, symmetrically to the left and to 
the right of the origin, since the rapid growth of the exponential term eventually prevails over the 
decay to zero of the inverse polynomial. The parameter v should only ensure that x*{j) € M. By 
relocating these minima x*{j) to itvr via rescaling, we obtain periodic functions gj G C°°[— 7r,7r), 
which decay in the time domain faster than any inverse polynomial. As a consequence, the high- 
pass filter Vj{x) = gj{x)v{x) is also quite smooth. The fact that some neighborhood B(0,5) is not 
contained in supp(i;) implies that the near-spikes of gj{x) at x = 0, a potential source of ripples 
in the time domain filter, disappear in the inverse Fourier transform of vj. Similar remarks apply 
to the low-pass filter Uj{x). The graph of the smoothed filter for j = 2 can be be seen in Figure 
5. The domain taken is [— vr, Sir) to illustrate that the periodic extension to M is quite smooth. 
Figure 4 displays the resulting time domain filter after numerical integration; moreover, for any 
fixed X, gj{2~^x) — )• g{x) as j — > oo. 

In the case of the fOU, the resulting filters gj, vj are displayed in Figure 6 for d = 0.25 and 
j = 2. The left and right plots illustrate the effect of multiplication by the high pass wavelet 
filter V. Note that, for d > 0, the filter gj^d, and thus gj, shows a singularity at the origin. This 
makes the role of the compact support high pass wavelet filter v quite important for the numerical 
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Figure 5: OU process: spectral filter 'gj for j = 2 {C, = 1). 



stability of the computation of the associated time domain filters. Analogously, the two plots in 
Figure 7 illustrate the effect of the low pass wavelet filter u. Even before multiplying by n, there is 
no singularity at x = 0, because the singularities in the individual terms gj, gj-i cancel out in the 
ratio gj{x)/gj-i{2x). The resulting low and high pass filters Uj, Vj in the time domain are shown 
in Figure 8. A generated sample path can be viewed in Figure 1. The remarkable persistence in 
the sample path, especially in comparison with that of the OU process, is due to the long range 
dependence of the fOU process when d > 0. 

As for the fGLE, smoothing is more challenging due to the presence of spikes in its spectral 
density (see Figure 2). Also, its more complicated functional form makes manipulation more 
difficult; however, we propose to continue to smooth via an exponential term due to the ensuing 
analytical simplicity. For large enough zbvr become approximate critical points of the component 
5j,7,d(3^)) thus dispensing with centering. In fact, from Table 2, the first order condition for the 
logarithm of the proposed filter in the range x > is that 

+ ^i^-'^' + = TW-'"' + Y722/3x^^-^ (2.15) 

The first and second terms on the left-hand side, and the first term on the right-hand side of 
(2.15) are close to zero for sufficiently large j. Thus, solving the resulting expression for x gives 
x*{j) = X* = vr as an approximation. Moreover, the approximation (1.8) holds as pointwise limit. 

Remark 2.2 The simplification behind (2.15) requires large enough j so that the critical point 
is close to vr. However, this depends on the values of the parameters ( and m. For instance, the 
high pass filter ui{x) for (Cni) = (2,1) is already quite smooth. However, for (C'm) = (10,1), 
ui{x) displays a kink at 7r/2, which is already much less visible in U2ix) (not shown). 

Remark 2.3 As discussed in Pipiras (2005), increasing the number of vanishing moments of the 
underlying MRA can improve the time domain decay of certain filters. Strictly speaking, the 
improvement in the decay depends on the specific Fourier domain form of the filter. (However, 
see Didier and Pipiras (2010), especially Remarks 6 and 7.) 

In the context of the present paper, we can give a slightly different explanation for the po- 
tentially positive effect of the number of vanishing moments. Increasing the latter has the effect 
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Figure 6: fOU process, high pass spectral filters at j = 2 = 1, d = 0.25). Left: gj. Right: vj 
(4 zero moments). 
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Figure 7: fOU process, low pass spectral filters at j = 2 {( = 1, d = 0.25). Left: gj{-)/gj^i{2-). 
Right: Uj (4 zero moments). 
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Figure 8: fOU process, time domain filters at j = 2 {( = 1, d = 0.25, 4 zero moments). Left: uj. 
Right: Vj. 



of increasing the regularity (Fourier domain smoothness) of the proposed filters. With a greater 
number of vanishing moments, Vj{x) becomes flatter over a wider vicinity of zero in the Fourier 
domain, which can more efficiently make up for a singularity or kink at the origin (e.g., for the 
fOU). In numerical studies, for = 4 or 8 vanishing moments and parameters values = 1, 
d = 0.25 we compared the time domain decay of fOU low- and high-pass filters uj and Vj. In 
general, for j = 2, 5, 8 and 11 the tail values (lags T = 31 through 40) of Uj and Vj under N = 8 
was on average of the order 10~^ below those obtained when = 4. Due to this small effect, we 
used = 4 all through the paper. 

Remark 2.4 All the time domain filters and covariance functions used in this paper were numer- 
ically calculated using the adaptive Lobatto quadrature method. For computational simplicity, 
we used Daubechies filters, instead of Meyer. In Matlab, the quadrature method is implemented 
via the quadl.m function. Section D provides a numerical study of the accuracy of the quadl.m 
function comparing the deviation of the numerically computed AR(1) and FARIMA filters from 
their closed form expressions. We also performed experiments with the adaptive Gauss-Kronrod 
quadrature, which is more suitable for functions with moderate singularities at endpoints. This 
method is implemented in Matlab via the quadgk.m function and is also studied in Section D. 
Note that, alternatively, the computation of Fourier transforms of functions displaying a singular- 
ity at the origin can be dealt with via a change of variables. (See Helgason et al. (2011), Section 
3.4.) 

3 Evaluating the simulation method 

In this section, we evaluate the accuracy of the simulation method. Most simulation techniques 
are supported by theorems that establish some sort of convergence, equality in law and so on. 
However, the finite sample performance can be disparate across methods in practice. One ap- 
proach is to use estimators to compare the simulation methods. Nevertheless, only the asymptotic 
distribution of estimators are available in most cases. The finite sample performance of estima- 
tors, e.g., bias, is then studied based on simulation, which creates a circularity. In view of this, 
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we study the performance of the methods relative to one another and compare with three other 
methods: simple iteration (OU process), Cholcsky and CME. 

Cholesky decompositions provide a classical and simple simulation method. If F is a target 
zero mean Gaussian stationary process with a given, known covariance matrix E over a finite set of 
time points N, a Cholesky decomposition of S = LL* is performed, where L is a lower triangular 

matrix, and vector Z of i.i.d. standard Gaussian variables is generated. Then V = LZ, as desired. 
Cholesky-based simulation is exact up to the accurate calculation of the covariance and can be 
implemented recursively (see Asmussen and Glynn (2000); see also Bardet et al. (2003), Craigmile 
(2005)). However, it is slow in terms of computational complexity: 0{N^), where N denote the 
length of the resulting stochastic vector. 

Another popular method is the CME (see Davies and Harte (1987), Wood and Chan (1994), 
Dietrich and Newsam (1997), Johnson (1994), Beran (1994), Asmussen and Glynn (2000), Percival 
and Constantine (2002), Craigmile (2003)). The algorithm involves embedding the covariance 
matrix in a non-negative definite circulant matrix of size M > 2(N — 1). This is computationally 
convenient, since the diagonalization of circulant matrices can be carried out by means of the Fast 
Fourier Transform (FFT), which has complexity 0(A^ log(A^)). Like Cholesky-based simulation, 
CME is exact. For a description of the CME, see Bardet et al. (2003), p. 582. 

Since the OU process can be simulated based on a simple loop, we choose this method to 
provide the baseline for the CME and the wavelet-based method. In the cases of the fOU and 
fGLE process, the baseline method is Cholesky-based simulation, since it is also a simple and 
exact procedure. A two-sample t statistic is used to assess the difference between the values of 
the estimator when generated by two of the methods. For the OU process, we evaluate the quality 
of the simulation based on the Yule- Walker estimator, whereas for the fOU and fGLE we use the 
Local Whittle estimation of the parameter d. The latter is of special interest in the framework of 
subdiffusion, since the Local Whittle is a good estimator for the subdiffusivity parameter a (see 
Didier et al. (2012)). 

For the OU process, the initial, exact step j = amounts to simulating through a simple 
loop an AR(1) process with parameter (p = and white noise variance ^— Based on 
both wavelet and exact simulation, we generated the OU process with parameters C = 1 
(7 = 1 over the interval [0, 2*^], with 2^^ points in each subinterval of length 1. Then, by sampling 
at the rate A = 2^^ (i-e., every 2^*^ points), the associated AR(1) process has parameter (/> = 
exp(-l • 2-3) « 0.8825, estimated by Yule- Walker over a time series of total length 2^^. In order 
to speed up the computations while preserving accuracy, the wavelet filters were truncated either 
at lag |r| = 40 or when a value below 10"^ is attained, whichever is first. In order to test the 
consistency of wavelet simulation for different values for the parameter C, and thus different filters, 
we also generated the OU process with parameters C = 2 and a = 1 over the interval [0, 2^] and 
sampled it at A = 2"^ and then with parameters C, = 1/2 and cr = 1 over the interval [0,2^] 
and sampled it at A = 2"^. Therefore, the associated AR(1) processes have the same parameter 
(f) = exp(— 2 • 2~^) = exp( — 1/2 • 2~^) ~ 0.8825. The simulation results, found in Table 3, suggest 
that the method is accurate when compared to iterative simulation and CME. The filters for 
C = 1/2 seem to be less accurate than those for C = 1,2. When the final scale is J = 4, the 
absolute value of the t statistic is above 4. However, when J is increased to 6 and 8, the latter 
drops below 2. This is indicative of increasing quality of the approximation as an increasing 
function of J. The same phenomenon is observed in Table 4, which displays results for C = 1 
when A = 2~^ and A = 2~^ (thus, cf) = exp(— 2~^) and ^ = exp(— 2~^), respectively). See also 
Table 5 additional results on the case when the time scries has total length 2^. 

For the fOU process the comparison is made over the integer time points 0,1,2, For 
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all simulations, filters were truncated at lag T = 40 giving entries on the order 10~^ at the 
point of truncation in the worst cases, typically in the low pass filters uj. We experimented 
with two different initializations: either via the convolution of the filter go := vq (i.e., at j = 0) 
with white noise or directly via CME with the associated autocovariance function calculated 
by means of numeric integration. The results can be seen in Table 6 for the parameter values 
d = 0.10,0.25,0.45. In practice, the disadvantage of initialization via direct convolution with 
white noise is that in principle it might require storage of fairly long filters when d > (i.e., under 
long range dependence). For this reason, go is truncated for some cases at T = 1, 200 in Table 6 
yielding a quite large overall length of the filter go at 2,401. The absolute value of the t statistics 
is less than 2 regardless of the initialization (convolution or CME), thus yielding results rather 
similar to CME for the simulation of the fOU. Though not displayed in the tables, simulation 
initialized with truncated filters at T = 400 or 600 seem to give rather similar results. 

For d < 0, the discrete time fOU process is not anti-persistent (see Corollary B.l). To 
evaluate the wavelet-based simulation procedure, we took the discrete time increment Y{n) = 
X{n) — X{n — 1) of the associated position process X{t) = Jq y(s)ds. The spectral density of 
Y{n) is, indeed, anti-persistent for d < (see Corollary B.l). For wavelet simulation purposes, 
the process X{-) was approximated by the simulated V{-) based on the expression (A. 21): for 
some fixed J, a sequence Vj^k was generated, and the approximation sequence 

L{2^-i)tJ , 

k=0 

was then calculated and sampled. The Cholesky sequences were simulated based on the covariance 
function of the process AX{t), 



r(2d + 2)sin(7r(d+l/2)) , 



EAX{s)AX{s + t) = ' / e 



1 



IX 



jdx. 



1^2 _j_ J.2 |-j;|2(i 



The results are also shown in Table 6 for different parameter values d = —0.10, —0.25, —0.45. 
Once again, similarly to CME, the absolute value of the t statistics is less than 2 in all cases. 

Table 7 displays a study of the accuracy of the wavelet-based simulation of fOU as a function of 
the finest scale J when d = 0.25, —0.25. Theoretically, as J — ^ oo, the quality of the approximation 
improves. However, the results show that the relative bias does not change much as a function 
of J for J = 2,4,6,8, 10. This potentially indicates that the quality of the simulation is already 
good enough at low values of J. 

For the fGLE, the filters displayed slower decay than those for the fOU process. For this 
reason, truncation was performed at lag |T[ = 80, which gave entries of the order 10~^ at the 
point of truncation in most cases, typically in the low pass filters, and entries of 10~^ only for 
the low pass filter ui for different values of d. For the same reason as for the fOU, the wavelet 
simulation was performed based on the expression (3.1) (see Corollary B.l). Also, the Cholesky 
sequences were simulated based on the covariance function of the process AX{t), 



EAX{s)AX{s + t) = c j e'*^ |- 



IX 



^ ' dx 

70 + 7i|2:|'^ +^2\x\^^ |xp<^ 



for an appropriate c > 0. The results are shown for different parameter values and final approxi- 
mation scales J in Table 8. In all cases, the absolute value of the t statistic is less than two and 
close to the corresponding value obtained from CME simulation. Also, in contrast with the OU 
process, increasing J does not seem to affect considerably the quality of the simulation. 
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Remark 3.1 Other simulation studies were carried out for the fOU, < d < 1/2, with filters 
calculated via the adaptive Gauss-Kronrod quadrature. The results were comparable to those 
obtained via adaptive Lobatto quadrature, so they are not shown. For all numerical integrals 
for all methods, the quadrature precision ranged from 10~^ to 10~^^. Whenever applicable, the 
length of the initial CME-simulated series for the wavelet method was 2^*^. 

4 Discussion 

Wavelet-based simulation methods have proven to be fast and efficient alternatives to FFT-based 
methods for rather large samples. They usually exhibit low computational complexity, since 
they are based on the Fast Wavelet Transform (see Percival and Waldcn (2000)). In this paper, 
we proposed an approximate wavelet-based simulation technique for two classes of continuous 
time anomalous diffusion models, the fGLE and the fOU. The proposed algorithm is an iterative 
method that provides approximate discretizations that converge quickly to the true sample path 
of the target process. The simulation technique involves an appropriate sequence of filters gj, 
j = 0,1,..., J, where is J the finest approximation scale chosen. The method then amounts to 
recursively generating an induced sequence of discrete approximations Vj, at scale j G N, where 
Vj = Qj * One can then show that 2'^/^V"j j^g^tj ~^ ■/ — > oo, in an appropriate sense, which 

naturally leads to an approximation of the position process X[t) via Riemann sums. As compared 
to previous works such as Didier and Pipiras (2008), this paper proposes a simulation procedure 
when the discretization of the target continuous time process at different scales does not have 
closed form in the time domain. Moreover, we propose smoothing procedures for the proposed 
filters as to speed their time domain decay, and thus minimize the border effect, which is always 
present in convolution-based procedures. 

While this method is approximate, it has several advantages such as: computational speed; 
approximations that converge uniformly over compact intervals almost surely; it is iterative; it is 
not intrinsically Gaussian. To study the performance of the wavelet-based simulation in compari- 
son to exact methods such as Cholesky and CME, we performed several Monte Carlo experiments. 
The simulation study measured the bias of well-established estimators when compared with real- 
izations from other methods. In most cases, the bias of the estimators when simulated using the 
wavelet method seems to lie within an insignificant distance from the exact methods. Therefore, 
the method is nearly as accurate with potentially reduced computational complexity compared 
with existing methods. 

A Proofs 

In this section, we discuss the adaptation of the original proofs in Didier and Pipiras (2008). 

We will need to replace Assumption 3 with the weaker Assumption 3', which allows for a kink 
in Gj and at ibTr, Assumption 4' replaces Assumption 4 for the sake of clarification. We 

remind the reader that the support of the Meyer scaling function (f){x) is contained in the interval 
[-47r/3,47r/3]. 

Assumption 3': 



Gj{x),Gj{x)-^ e C[-47r/3,47r/3] nC2([-47r/3,47r/3]\{±7r}) 



(A.l) 



max sup sup \Gj{xy\ < oo 

'=-l,lj>jo \x\<4:Tt/3 



(A.2) 
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max max sup sup 

p=-l,l k=l,2 \x\jLn,\x\<'ln/3 



max max sup lim 

p=-l,l fe=l,2 j>jg a;->-7r+,— TT- 



gk 



Qk 



< oo 



< oo 



(A.3) 
(A.4) 



Assumption 4': ^ is twice differentiable in M\{0} and, for large |x|, 



X 



< 



const 



|fc+i ' 



A; = 0,1, 2. 



(A.5) 



Assumption 4' is clearly satisfied by the spectral filters g{x) of all the processes considered in this 
paper, so we turn to Assumption 3'. 

In Lemma A. 2, we show that the conclusion of Lemma 1 in Didier and Pipiras (2008) still 
holds if Assumption 3 is replaced with Assumption 3'. For the reader's convenience, we discuss 
the modification of the proof. Before that, we establish a simple auxiliary lemma. 

Lemma A.l Let f,g be two real functions in C^{a,b) such that f, g' are bounded over [a,b] and 
the limits ^^x-^a+,b- f{x)9{x) exist. Then the formula for integration by parts holds, i.e., 

lim fib - eb)g{b - £&) - lim /(a + ea)g{a + £„) = / [f\x)g{x) + f{x)g' {x)]dx. 

£6->0+ £a-*-0+ Ja 

Proof: Since f,g & C^{a,b), then for small ea,e6 > 0, 

rl>-£b 

f{b - eb)g{b - £b) - f{a + £a)gia + £„) = / [f'{x)g{x) + f{x)g' {x)]dx. 

■J a+Ca 

By the remaining assumptions and taking the limits lim£^_^o+) li™£6->0+) the claim follows. □ 



Lemma A. 2 Under Assumptions 3' and 4' , 

|2-^/2$j(2-%)|, \2-^/^^^{2-^u)\ < Y^q^, ueR (A.6) 

C2-J'/2 

\^'{2-^u)\<-—^, ueR, (A.7) 

where $j is as in (1.11) , and , are defined in the Fourier domain as $-^(x) = 

Gj{2-ix)-^2-^l'^4){2-^x), ^^{x) =g{x)2-3/'^i){2-^x). 

Proof: We first look at . Consider initially x > 0. From the properties of the Meyer MRA, 
Gj{x)~^(l){x) also satisfies Assumption 3'. Thus, by Lemma A.l, 



47r/3 47r/3 piux piTr/Z tux a 

e''^''Gj(x)-^(j)(x)dx= GAx)-^4>{x)- — {GAx)-^(t>ix))dx 

■K lU ./_ lU OX 



/ ^ ^3 

J-rv 



-i (e»-G,(7r)-^^(7r) + J^' e^^^ -^{G,{x)-'^{x))dx) . (A.8) 
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In turn, again by Lemma A.l, the integral on the right-hand side of (A.8) is 



TT iu dx 



/•47r/cS tux fiZ 



c\ /*A.TT 1 3 tux 



As for X < 0, an analogous expression holds for 



/ ^ e'"^Gj(x)-^0(x)dx. (A.IO) 

i-47r/3 

On the other hand, 

J — TT 

= ^(e^""G,(7r)-^^(7r) - e— -G,(-7r)-^^(-;r) - J^^^ ^{G,{x)-'^{x))dx) . (A.ll) 
Once again by Lemma A.l, the integral on the right-hand side of (A.ll) is 

li^'^\'}^.^x^''^^^^-"^^^^^ 

~ j\ ^''''^iGj{x)-'^{x))dx) . (A.12) 
Consequently, by adding together (A.8), (A. 9), (A.IO), (A.ll), (A.12) and by Assumption 3', 

r47r/3 



2-j/2$j(2-J-^)| = ( / +/ +/ )e^""G,(x)-i0(x)dx 

2vr 7-47r/3^ 



< 



27rit2 



for a constant C that docs depend on j, which gives the inequalities in (A. 6) for The remaining 
inequality, for ^j, can be obtained by a similar procedure. 
To show (A. 7), start from 

2i/2 



^^{2~^u) = — [ e''''^g{2^x)^{x)dx, « G M. 
27r 



Since the only possible singularity of 'g is at the origin by Assumption 4' and suppjV'} C {27r/3 < 
|x| < Stt/S}, the same argument as in the proof of Lemma 1 in Didier and Pipiras (2008) applies, 
thus yielding (A. 7). □ 

The following lemma is a claim in the proof of Proposition 1 in Didier and Pipiras (2008). 
Since we use filters under slightly different assumptions in this paper, we state it and provide a 
more detailed proof. 

Lemma A.3 Let u, v be Meyer CMFs such that u{x) G C^[— 7r,7r). Then, under Assumptions 3' 
and 4' , 

Kkl\Hk\<o(^j^y (A.13) 
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Proof: The Meyer low-pass filter u satisfies suppjii} C {|a;| < 2it/3}. On the other hand, the 
possible kinks of Gj^i{x) in [—An/ 3, in/ 3] lie at ztyr by Assumption 3'. Therefore, the poten- 
tial kinks of Gj{2x) in [— 27r/3, 27r/3] lie at ±7r/2. However, those points lie outside suppjS}. 
Therefore, we can write 

^V-27r/3 J-tt/2 Jn/2 ' Gj[2x) 

and again by Assumption 3', one can use the same type of argument as in the proof of Lemma 
A. 2 to establish (A. 13) for uj^k- 

As for Vj^k, dj+i{x) is smooth except possibly at the origin by Assumption 4'. Since Vj{x) = 
gj+i{x)v{x), supp{w} C {tt/S < \x\ < Btt/S} and the fact that v{x) G C^[— 7r,7r), then by applying 
integration by parts twice we arrive at the claim. □ 

We now describe the necessary modifications to the remaining claims in Didier and Pipiras 
(2008): 

• Theorem 2, section 5, expression (5.13) we still have that 

Fm{x) := ( me-''""') .4ff-2-^/2^(2-^x) ^ 9^ix), m ^ oo 

k=-m ^J^'^ ^> 

in L^(M), since 5j,jte~''^^ converges to ^j(x) in tt, tt) and, by Assumption 3', 

g,^(^Jx) bounded on the compact support of (f){2~'^x). 

• Proposition 1, section 6, expression (6.5): as an immediate consequence of Lemma A. 3, 

Uj{x) = Gj+i{x){Gj{2x))-^u{x) e L^[-7r,7r), Vj{x) = Gj+i{x)g{2^+'^x)v{x) G L'^[-tt, it). 

(A.14) 

By a similar proof as that for Lemma A. 3, we also conclude that 

n^(x) = Gj+i{x)-^Gj{2x)u{x) G L'^[-Tr,Tr), v^{x) = Gj+iix)-^g{2j+^x)-^v{x) G L'^[-tt,tt) 

(A.15) 

We now show that the proposed filters satisfy Assumption 3'. Let 

fj{x) = Gj,^,d{xf. (A.16) 

Then fj{x) satisfies (A. 2) (i.e., fj and its inverse is uniformly bounded over |x| < 47r/3 and large 
j). Moreover, since 

G'j,^,a{^) = lMx)-'/'f'j{x), G;^,,(x) = i(- i/,.(a;)-3/2/;.(x)2 + /,.(x)-V2/;.'(a;)), 
then it suffices to look at fj{x) and f'j (x). 

Lemma A. 4 Let Gj^(^{x), Gj^dix) and Gj^^^d{x) he as in (2.12), (2.13) and (2.14), respectively. 
Then all these filters satisfy Assumptions 2, 3', and 5. 
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Proof: It is clear that Assumptions 2 and 5 are satisfied in all cases, so we focus on Assumption 
3'. 

The argument is straightforward for Gj^(^(x) and Gj^d{x)- As for Gj^^^di^), let fj{x) be as in 
(A. 16). Then it suffices to show that it satisfies (A. 3) and (A. 4) (with /. in place of G.). Note that 
fj{x) = l,\x\<Tr, and thus the first and second derivatives are trivial in this range. Without 
loss of generality, we now only look at the range tt < x < 47r/3. In this case, 

^ 70 + 7i2J/3(27r - x)^ + 7222i/3(2;r - x)2/3 ' 

Thus, 

f'jix) = {^i2^^Px^-^ + ^22^^^2Px^^-^){jo + 7i2^'^(27r - x)^ + 7222^'^(27r - xf^y^ 

+(-1)2(70 + 7i2^'3x'' + ^22^^^x^^){jo + 7i2^'^(27r - xf + 7222^'3(27r - xf^)-^ 

•(7i2J'''^(27r - x)'^-^ + 7222^'^2^(27r - xf>^-^) 

=: aj{x)bj{x) + Cj{x)dj{x)ej{x). (A.17) 

As for the first term in the sum (A.17), note that, for any e > 0, for large j, \wj\ := |7o/22-'^ _|_ 
7i/2-^'^(27r — x)^| < e uniformly in x over tt < x < 47r/3. Moreover, since 72 = m? > 0, then 
72(27r — x)2^ attains its (constrained) minimum at x = 47r/3. Therefore, 

\bj{x)\ < |72(27r - 47r/3)2^ - \wj\\-^2-^^^, 

whereas 

|a,(x)| < 22^/^(^^7r/^-i + |72|2/37r2/3-^). (A.18) 

Now consider the absolute value of the second term in the sum (A.17). By a similar reasoning, 

|c,(x)d,(x)e,(x)| < 22./^ + t^(|)^^,(|[)^)2(-2)^^•/^|72(2. - 4./3)^^ - 

•22^-/3 (^/3(27r - 47r/3)^-i + 722/3(27r - An/Sf^-') . (A.19) 

By (A.18) and (A.19), |/j(x)| (and thus also {fjix)~^)') is uniformly bounded over |x| < 47r/3 
and large j, i.e., it satisfies (A.2). 

We now turn to the second derivative. We obtain 

/;(x) = [{ji2^''Pif3 - l)x'^-2 +7222^-/32/3(2/3 - l)x'^-^)bj{x) 

+aj(x)(-l)2(7o+7i2^''3(27r-x)'^+7222^'^(27r-x)2/^)-2.(7i2^'^;3(27r-x)^-^+7222j'/32^(27r-x)2'3-i)" 

+ 

(7i2J''5^x^-i + 7222^-^2^x2/3-1) [d,(x)ej(x)] + c,-(x)[d,-(x)e,-(x)]', 

where 

[dj{x)ej{x)]' = (-2)(7o + 7i2^/3(27r - xf + 7222^'3(27r - xf^^)-^ 

(7i2^/3/3(27r - x)/'-i(-l) + 7222j/52/3(27r - x)2/3-i(-l))ej(x) 

+d,-(x)(7i2^'^/3(;9 - l)(27r - xf-\-l) + ^22^^Hl3{2l3 - \){2it - xf^-^{-l)). 

By a similar reasoning to that for f'jix), we therefore conclude that \ fj{x)\ (and thus {fj{x)~^)") 
is uniformly bounded over |x| < 47r/3 and large j. □ 



21 



Proposition A.l Let {V{t)}t>o be the velocity process for the fGLE (1.2), (1.4) or the fOU 
(1.5). Then, 

sup\2-^/^Vj,2Jt] - y{t)\ < Ai2--^'' a.s. (A.20) 
teK 

for some v E {0,1), where K is compact interval and Ai is a random variable that does not depend 
on J. 

Proof: By the proof of Proposition 2, section 7, in Didier and Pipiras (2008), we have to show 
that V satisfies Assumptions 2, 5, and a Holder condition, and make use of Assumption 3' instead 
of Assumption 3. In fact. Assumptions 2 and 5 are satisfied by Lemma A. 4. In Lemma B.l we 
estabhsh the Holder condition (B.2) for some i/ G (0, 1) in the cases of the fGLE and fOU. Finally, 
we can bound |2~"^/^$j(2~'^i;)| based on Lemma A.2, which is a consequence of Assumptions 3' 
and 4'. The latter is satisfied for the spectral filters g{x) of the processes in question, whereas the 
former also holds in view of Lemma A. 4. Thus, the claim follows. □ 



Remark A.l Note that, in the case of this paper, Assumptions 6, 3*, 5*, 7 of Didier and Pipiras 
(2008) are not used since the Holder continuity order given by Lemma B.l is v E (0, 1). 

The next proposition shows that the position process X can also be uniformly approximated 
by a partial sum process defined based on the discrete approximation to V. Without loss of 
generality, we show it for T G N. 

Proposition A.2 Let {V{t)}t>o be the velocity process associated with the fGLE (1.2), (1.4) or 
the fOU (1.5), and let T G N. Consider the sequence {Vj^k}k=o ...2Jt-1j where each term can be 
interpreted as Vj^k = Vj |^2jJ^j • Then 



sup 
te[o,T] 



X{t)- Yl 2^/2K7,fc^| <A2-^^ ^^£(0,1), 

k=0 



(A.21) 



where the random variable A does not depend on J. 



Proof: By Lemma B.l, X{t) is well-defined as the integral (B.l). Fix J > and form an 

associated partition < ttt^ > of [0,rl. Assume that T > t > 1. The case where t < 1 

) k=0,...,2JT-l ^ ' ~ ~ 

can be handled similarly. 



On the one hand, 

. L{2-'-i)tJ 



L(2^-l)tJ 



' ~ k=0 k=0 

J, L(2^ - m + 1 



fe=0 



1 



< Ai2- 



2J 



(A.22) 



where the last inequality follows by Proposition A.l. 



On the other hand, note that, since t > 1, then 2-'^*"^"'"'^ — Therefore 



V 



k=0 



k/2J 



(f) 



ds 



22 



L(2^-l)tJ (fc+i)/2J 



L{2-'-i)tJ 



r-{fc+l)/2-' 



ds 



Ik/2-' 



< max sup 

fc=0,l,...,L(2-^-l)tJ 

< max sup A2 

fc=o,i,...,L(2-'-i)tJ sgL4,^J 



L(2^ - m + 1 



2-^ 



s — 



2-^ 



[(2-^ - l)tj + 1 
2-^ 



< ^32-^'^ 



(A.23) 



by the Holder condition (B.2), where A2, are random variables that depend only on T. Also, 
by the Holder continuity of V (Lemma B.l), 



f V(s)ds < sup \V(s)\(t 

Jl(2J-2W±l ~ se[o,T] ^ 



[(2-^ - 1)^J + 1 
2^ 



) < sup \V{s)\ C2--^ (A.24) 



for some constant C > 0. 

Now taJce sup(g[o,T] both sides of (A.22), (A.23), (A.24). The claim follows from the triangle 
inequality. □ 



B Auxiliary results 

In this section, wc develop the Fourier domain integral representations for X. We first establish 
that X can be regarded as the integral of V . 

Lemma B.l Let {V{t)}t>Q he the velocity process for the fGLE (1.2), (1.4) or the fOU (1.5). 

Let ^ 

X{t) := I V{s)ds, (B.l) 

where the integral (B.l) is taken in the Lebesgue sense. Then (B.l) is well-defined in the sense 
that there exists a process r){t) which is equivalent to V{t), and which satisfies the Holder condition 

\v{t) - r]{s)\ < A\t - si" a.s. (B.2) 

for some random variable A that only depends on [0, T] and some v G (0, 1). 
Let {X(t)}t>o be the position process associated with {V(t)}t>o. Then 

X{t) = (^-^)g{x)B{dx), 0<d<-, (B.3) 
\ IX / z 

holds, where g{x) is a spectral filter and B{dx) is given in (2.1). 

Proof: We only look at the fGLE, since the argument for the fOU can be developed along the 
same lines. 

We first show (B.l). From Cramer and Leadbetter (1967), pp. 181-182, the conclusion follows 
from verifying the condition 

/•oo 

/ x^''\og{\ + x)\g{x)\^dx <oo 
Jo 
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from some v € (0,1). Let £ > 0. Since |5(-)P in (2-5) is continuous, the only potentially 
problematic points arc the origin or oo. As x ^ 0"*", |^(x)p ^ jxp'^, so x^'^ \og{l + x)\g{x)\'^dx < 
oo for > 0. As X ^ oo, |5(x)|2 ~ x~'^^^~'^\ so /^°° x'^^ log(l + x)\g{x)\^dx < oo for z/ < + 1/2. 

To show (B.3), note that s,s' > 0, E\V{s)V{s')\ < ^E{V{s)y^E{V{s')y, which is finite 
and constant, by stationarity. Now apply Fubini's Theorem and formula (B.l). □ 



The following is a corollary to Lemma B.l. 



Corollary B.l Let {X{t)}t>o be the position process associated with the fGLE (1.2), (1.4). De- 
note the discrete-time first difference process by := AX(n) = X{n + 1) — X{ri), n G N U {0}. 
Thus, its spectral density fy is 



frix) = c{df 



( 




2 




ix 





^/«(x)|x| — imx 



-2d 



kei\{o} 



Ck{x + 2Trk)\x + 27rA;|-2'^ - im{x + 2Trk) 



\x + 



27rfc|-2('^+i)). 



Let {y{t)}t>Q o-nd {X{t)}t>Q be the velocity and position processes, respectively, associated with 
the fOU (1.5). Then the spectral density of the discrete time process {T^(t)}t=o,i,... is 



^ . , r(2d + 2)sin(7r(d+l/2)) / 1 1 ^ 



kez\{o} 



C^ + {x + 2A;7r)2 \x + 2A:7r| 



(B.4) 

Denote the discrete-time first difference process by := AX(n) = X(n + 1) — X(n), n G Nu{0}. 
Thus, its spectral density fy is 



frix) = 



r(2d + 2) sin(7r(d+l/2)) 



27r 



1 - e" 



IX 



1 



1 



+11 -e 



-ixi2 



ke7\{0} 



C2 + (x + 2/c7r)2 \x + 2/c7r|2'^+2 / ' 



C2+x2 \x\'^d 

X G [— TT, vr). 
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C Tables 



Table 3: Yule- Walker estimation of the AR(1) parameter ^ = exp(— 1 • 2~'^) = exp(— 2 • 2""') = 
exp(— 1/2-2~^) 0.8825 (wavelet filters cut off at lag |r| = 40 or at value 10~^, time series length 
2^^ for all methods) 



(p w 0.8825 




4> 


s 




\t\ statistic 


wavelet (C = 1, J = 


= 6) 


0.88016217 


0.01047668 


5000 


0.00300676 


wavelet {( = 2, J = 


= 8) 


0.88017601 


0.01072520 


5000 


0.27145002 


wavelet (C = 1/2, J 


= 4) 


0.87919411 


0.01072708 


5000 


4.34765925 


wavelet (C = 1/2, J 


= 6) 


0.87978815 


0.01064935 


5000 


1.55891496 


wavelet (C = 1/2, J 


= 8) 


0.88023961 


0.01058475 


5000 


0.57449953 


CME 




0.88051697 


0.01055337 


5000 


0.60412089 


iterative 




0.88011831 


0.01052941 


5000 





Table 4: Yule- Walker estimation of the AR(1) parameter cp (wavelet filters cut off at lag \T\ =40 
or at value 10~^, time series length 2^^ for all methods) 



4> 0.7788 






s 




\t\ statistic 


wavelet (^ = 1, J = 


6) 


0.77618232 


0.01385836 


5000 


2.50237872 


wavelet (C = Ij J = 


8) 


0.77670932 


0.01388566 


5000 


0.61614593 


iterative 




0.77688169 


0.01408879 


5000 




(f) ^ 0.9394 






s 


A^ 


\t\ statistic 


wavelet (C = 1; J = 


6) 


0.93651069 


0.00790679 


5000 


3.81308606 


wavelet (C = 1) J = 


8) 


0.93682813 


0.00784243 


5000 


1.80796495 


iterative 




0.93711216 


0.00790679 


5000 





Table 5: Yule- Walker estimation of the AR(1) parameter (j) (wavelet filters cut off at lag |T| = 40 
or at value 10~^, time series length 2^ for all methods) 

(j) « 0.8825 ^ ~s N |t| statistic 

wavelet (( = 1, J = 6) 0.87357855 0.02195045 5000 0.64730215 

wavelet (( = 2, J = 8) 0.87377009 0.02183448 5000 0.20701114 

wavelet (( = 1/2, J = 4) 0.87288940 0.02214156 5000 2.21835173 

iterative 0.87386073 0.02164170 5000 
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Tabic 6: fOU: Local Whittle estimation of o! (C = 1) (wa.vclct filters cut off at lag |r| = 40, time 
series leiiglli 2" for all methods) 



d = 0.10 






d 


s 


N 


\t\ statistic 


wavelet (qq length 2 x 400 + 1, 


J = 


6) 


0.12106833 


0.09259765 


5000 


0.01292907 


wavelet (CME at j = 0, J = 


= 6) 




0.11820450 


0.09406420 


5000 


1.65850145 








nil nnoi f\A 
(J.liyyzl(J4 




5000 


{j.74:24:2766 


Cholesky 






0.12130840 


0.09308411 


5000 


~ 


d = 0.25 








s 


N 


\t\ statistic 


wavelet {go length 2 x 1200 + 1, 


J = 


:6) 


0.27379919 


0.09391947 


5000 


1.19695166 


wavelet (CME at j = 0, J = 


= 6) 




0.27623054 


0.09390228 


5000 


0.10736144 


CME 






0.27474721 


0.09585494 


5000 


0.68122236 


Cholesky 






0.z7d03043 


0. 09248460 


5000 




d = 0.45 






d 


s 


N 


\f.\ statistic 


wavelet (2 x 1200 + 1, J = 


6) 




0.48484741 


0.09462295 


5000 


1.04322102 


wavelet (2 x 1400 + 1, J = 


6) 




0.48596691 


0.09366874 


5000 


1.64788134 


wavelet (CME at j = U, J = 


-- 6) 




U.4ozoo4z4 


U.U9478310 


5000 


U. 26842086 


CME 






0.48504209 


0.09394773 


5000 


1.15107347 


1 1 

Cholesky 






0.48288866 


0.09313095 


5000 




d = -0.10 








s 


N 


\t\ statistic 


wavelet (CME at j = 0, J = 


= 6) 




— U.Ud958o97 


0.09376265 


5000 


1.87458648 


CME 






—0.07254477 


0.09363777 


5000 


0.28587914 


Cholesky 






-0.07307713 


0.09257674 


5000 




d = -0.25 






d 


s 


N 


|i| statistic 


wavelet (CME at j = 0, J = 


= 6) 




-0.21858996 


0.09392950 


5000 


0.93216739 


CME 






-0.21875428 


0.09171082 


5000 


1.03200374 


Cholesky 






-0.21684272 


0.09350826 


5000 




d = -0.45 






d 


s 


N 


\t\ statistic 


wavelet (CME at j = 0, J = 


= 6) 




-0.40363919 


0.09449238 


5000 


1.07109047 


CME 






-0.40688947 


0.09378554 


5000 


0.64835641 


Cholesky 






-0.40566672 


0.09480350 


5000 
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Table 7: fOU: Local Whittle estimation of c? ((" = 1), comparison across values of J (wavelet filters 
cut off at lag 1^1= -10, lime series length 2" for all methods) 



d = 0.25 






d 


s 


N 


\t\ statistic 


wavelet (CME at j = 


0, J 


= 2) 


0.27497237 


0.09367323 


5000 


0.56835559 


wavelet (CME at j = 


0, J 


= 4) 


0.27473605 


0.09096327 


5000 


0.70556038 


wavelet (CME at j = 


0, J 


= 6) 


0.27430310 


0.09215599 


5000 


0.93550788 


wavelet (CME at j = 


0, J 


= 8) 


0.27408226 


0.09386958 


5000 


1.04538369 


wavelet (CME at j = 


0, J -- 


= 10) 


0.27750365 


0.09482862 


5000 


0.78643924 


Cholesky 






0.27603043 


0.09248460 


5000 




d= -0.25 






d 


s 


N 


\t\ statistic 


wavelet (CME at j = 


0, J 


= 2) 


-0.21832655 


0.09358734 


5000 


0.79308646 


wavelet (CME at j = 


0, J 


= 4) 


-0.21692926 


0.09455304 


5000 


0.04601620 


wavelet (CME at j = 


0, J 


= 6) 


-0.21911564 


0.09132082 


5000 


1.22965544 


wavelet (CME at j = 


0, J 


= 8) 


-0.21530451 


0.09444787 


5000 


0.81837755 


wavelet (CME at j = 


0, J -- 


= 10) 


-0.21680078 


0.09295220 


5000 


0.02249260 


Cholesky 






-0.21684272 


0.09350826 


5000 
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Table 8: fGLE: Local Whittle estimation of d {( = 2. m = 1), comparison across values of J 
(wavelet filters for AX{n) ciil olf al lag \T\ = 80: time series of length 2" for all methods) 



d = 0.10 



d 



N 



\t\ statistic 



wavelet (CME at j = 0, J = 6) 0.12044083 0.09245007 5000 1.05431590 

wavelet (CME at j = 0, J = 8) 0.12241433 0.09102200 5000 0.02227624 

wavelet (CME at j = 0, J = 10) 0.11932462 0.09357999 5000 1.65282860 

CME 0.12221790 0.09220211 5000 0.08515547 

Cholesky 0.12237381 0.09088329 5000 



d = 0.25 



d 



N 



\t\ statistic 



wavelet (CME at j = 0, J = 6) 0.27857299 0.09217255 5000 0.25582876 

wavelet (CME at j = 0, J = 8) 0.27998323 0.09268102 5000 0.50777309 

wavelet (CME at j = 0, J = 10) 0.28029380 0.09627634 5000 0.66273074 

CME 0.27818405 0.09112314 5000 0.46947509 

Cholesky 0.27904459 0.09144735 5000 



d = 0.45 



d 



N 



\t\ statistic 



wavelet (CME at j = 0, J = 6) 0.43156246 0.10538368 5000 0.78943448 

wavelet (CME at j = 0, J = 8) 0.43199573 0.10660877 5000 0.58088823 

wavelet (CME at j = 0, J = 10) 0.43166818 0.10678505 5000 0.73448769 

CME 0.43338219 0.10403330 5000 0.07274950 

Cholesky 0.43322954 0.10579033 5000 



D A study of the accuracy of numerical integration 

In Table 9, we study the accuracy of quadl.m by comparing it to the closed-form 
V-i = ^''"(1 - <l>^~'"T^dx = <^lo>o}, -K < 1. 

For |(/)| = 0.5, we choose T = 35 because of the early decay of the filter to values comparable to 
machine precision. 

Table 9: Numerical accuracy, AR(1) filters: quadl.m versus closed-form (cut off at lag |T|, 
deviation^ = ^j*"^^^ — V't) 





|r| 


mean abs. dev. {\t\ < \T\) 


-0.95 


100 


6.80 X lO"^'^ 


-0.50 


35 


8.33 X 10-^'^ 


0.50 


35 


3.04 X 10-^'^ 


0.90 


100 


3.22 X 10"^'^ 


0.95 


100 


3.73 X 10"^'^ 



In Table 10, we study the accuracy of quadl .m by comparing it to the closed-form 

The value p gives the radius of the ball around the singularity. The results show a high degree 
of accuracy in all depicted instances of AR(1) and FARIMA(0,d,0). As expected, the integration 
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Table 10: Numerical accuracy, FARIMA filters: quadl.m versus closed- form cut off at lag |r| 
deviation^ = -0^"^^^ — ipt) 

d p \T\ mean abs. dev. {\t\ < \T\) 



-0.45 10-"^^ Too 2.26 x 10"^* 

-0.35 10-^2 100 2.39 x 10"^^ 

-0.25 10-12 3.01 X 10-1^ 

-0.10 10-12 100 4.24 X 10-1^ 

0.10 10-12 100 5.56 X 10-12 

0.25 10-12 100 3.92 x 10-i° 

0.35 10-12 100 6.62 x 10-^ 

0.40 10-12 100 2.71 X 10-^ 

0.45 10-11 100 3.92 x 10-^ 



procedure becomes less accurate in the case of FARIMA as d approaches 0.5. For d = 0.45, for 
instance, we picked p = lO-n due to convergence problems. To compute the covariance function 
of FARIMA(0,d,0) 

I ^^-r(d)r(i-d)r(i-d + i)^{^-^°>' 2'2J\-f°^' 

quadl.m required a wider radius around the singularity at zero as d approaches 0.5 (not shown). 
The command quadgk.m is suitable when the singularity at zero behaves like for p > —0.5. 
The results are displayed in Table 11. The deviation for the lag f = is shown separately because 
it is the only one of a different order of magnitude. 



Table 11: Numerical accuracy, FARIMA autocovariance: quadgk.m versus closed-form (cut off at 

lags 1 and T, deviation^ = 'y^^^'l^ — It) 

d T mean abs. dev. (1 £ ^ < T) abs. dev. (t = 0) 

0.10 100 9.66 X 10-12 1.95 x 10^2 

0.15 100 1.98 X 10-11 4.88 x 10-2 

0.20 100 2.02 X 10-1° 9.87 x 10-2 

0.25 100 1.36 X 10-13 1.80 x 10"! 
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E Pseudocode 



Wavelet-based simulation over the interval [0,1] 



Input: the parameter values of the stochastic process, the final approximation scale J, 
the number of zero moments of the underlying wavelet basis; 

Step 1: numerically calculate the Fourier transform of Uj{x), Vj{x), j = 0, 1, . . . , J, to obtain 
the time domain filters Uj , Vj . For simplicity, let L be the constant length of Uj , Vj ; 

Step 2: numerically calculate the autocovariance function ro of the discrete time process 
induced by the filter go{x) at j = 0; 

Step 3: generate an initialization sequence Vo of size L + 1 via CME with the autocovariance ro; 

Step 4: generate a sequence of standard Gaussian white noise Sj of the same length 

(2^ — 1) + (L + 1) of Vj. Insert zeroes between the entries of both vectors Vj, £j. Convolve the 

resulting vectors with the low- and high-pass filters Uj, Vj to obtain 

Uj* t2 Vj, Vj* t2 Ej, respectively. 

Step 5: drop the 2{L — 1) entries of the vectors Uj* t2 Vj, Vj* t2 £j which are affected 
by the boundary effect (filter truncation) and add together the two resulting vectors 
to obtain the vector Vj+i of length {2^+^ - 1) + (L + 1). 

Step 6: stop if j = J. The resulting simulation vector Vj has size (2-^ — 1) + (L + 1) > 2^. 
Step 7 : else set j j + 1 and goto Step 4. 



Remark E.l Simulation over the time interval [0, 2^], K ^N, can be performed in an analogous 
fashion. 
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